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ABSTRACT 

Observations show a clear vertical metallicity gradient in the Galactic bulge, which is often taken 
as a signature of dissipative processes in the formation of a classical bulge. Various evidence shows, 
however, that the Milky Way is a barred galaxy with a boxy bulge representing the inner three- 
dimensional part of the bar. Here we show with a secular evolution N-body model that a boxy 
bulge formed through bar and buckling instabilities can show vertical metallicity gradients similar to 
the observed gradient, if the initial axisymmetric disk had a comparable radial metallicity gradient. 
In this framework the range of metallicities in bulge fields constrains the chemical structure of the 
Galactic disk at early times, before bar formation. Our secular evolution model was previously shown 
to reproduce inner Galaxy star counts and we show here that it also has cylindrical rotation. We use 
it to predict a full mean metallicity map across the Galactic bulge from a simple metallicity model 
for the initial disk. This map shows a general outward gradient on the sky as well as longitudinal 
perspective asymmetries. We also briefiy comment on interpreting metallicity gradient observations 
in external boxy bulges. 

Subject headings: Galaxy: structure — Galaxy: bulge — Galaxy: abundances — Galaxy: kinematics 
and dynamics — galaxies: evolution — methods: numerical 



1. INTRODUCTION 

About half of edge-on disk g alaxies contain box y 
or peanut-shaped bulges (BPBs, iLiitticke et al.l [2OOOI ). 
Their photometric and kinematic properties are con- 
sistent with predictions from disk galaxy simulations 
in which a BPB formed through bar and buckling 
instab ilities (iBureau fc Freemaiil 119991 : iDebattista et al.l 
120041: IBureau fc Athanassoulal 120051 ). The Galactic 
bulge shows many characteristics of a BPB formed in 
this wa y: a boxy shap e in projection (|Dwek et al.lll995l : 
[Skrutskic et al. 20 06[): a triaxial density distribution 
(p3inncy e t al. 1997t~ fL6pez-Corredoira et a l.l ^2005) with 
a den se inner, rounder component ([Gonzale z "eFall 
I2012t iGerhard &: Martinez- Valpuesta ' ^2012, he reafter 
GMyi2); an X-shaped s t ructure jM cWiUiam & Zoccah 
2010[ iNataf et al.l 120101: iNess et aL 20 12at iLi fc Shen 
2OI2I): cylindr ical rotation ( Howard et al.l 120091: 
Shen et all I2010D: and a transition to an o uter pla- 



nar bar (jMartinez- Valpuesta fc Gerhardl[20Tll hereafter 
MVGll), seen in star count observations as 'lon g bar' 
(|Beniamin et al.ll2005l: iCabrera-Lavers et al.ll2007D . 

The Galactic bulg e is also known to consist of pre- 
dominantly old stars ([Zoccali et al.l l2003': Cla rkson et al.l 
[2OO81) with a broad, asymmetric metallicity distribu- 
tion function (MDF ) (IRichlll990t llbata fc Gilmord [19951: 
IZoccali et al.l I2003[) . The data clearly show a vertical 
metallicity gradient, such that the more metal-rich part 
of the MDF thins out towards hig h latitudes (b < —4°, 
Minniti et all 119951: iZoccah et a"Lr[200 8: Go nzalez et alj 
20111) . At latitudes lower than Baade' s window (—4° < 
b < 0°) no vertical gradient is seen (IRich et al.l [2007L 
I2OI2D . These properties are not obviously reconciled with 
a disk origin of the bulge. Because bars tend to mix stars 
from different radii and so wash out preexisting metal- 
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licity gradients (jFriedli et al.lll99"^ , the observed metal- 
licity gradient in the Galactic bulge has long been taken 
as a s ignature for a classical bulg e in the Milky Way. In- 
deed, IBekki fc Tsuiimotol (|2011[ ) found that they could 
not reproduce the observed gradient with a simulated 
boxy bulge which had evolved from a pure disk with an 
initial metallicity gradient. In order to obtain the lower 
mean metallicities at high 6, they needed to include an 
additional metal-poor thick disk in their initial model. 

Here we show with the help of a suitable boxy 
bulge/bar simulation, that contrary to widespread be- 
lief boxy bulges can have metallicity gradients similar to 
those observed in our Galaxy and in other boxy bulges, if 
they evolve from disks with similarly steep radial metal- 
licity gradients. 

2. N-BODY MODEL FOR THE GALACTIC BAR AND BOXY 

BULGE 

The simulation used in this work is that already an- 
alyzed in MVGll an d GMV12, with similar ch aractcr- 
istics as described in [Martinez- Valpuesta et al.l p006,) . 
The model evolved from an exponential disk, with Q = 
1.5, initial radial scale length of /i = 1.29 kpc and verti- 
cal scale height of hz = 0.225 kpc, embedded in a dark 
matter halo. Following bar and buckling instabilities it 
developed a prominent boxy bulge which is already re- 
laxed at ^ 1.5 Gyr. As in GMV12 we consider the sim- 
ulated galaxy at time Tb ~ 1.9 Gyr when the bar has 
resumed its growth through further angular momentum 
transfer to the halo, and we vertically symmetrize the 
particle distribution to optimize the resolution. 

We showed previously that this model reproduces dif- 
ferent Milky Way star count observations, both near the 
Galactic plane for the long bar out to / 27° (MVGll), 
and in the inner boxy bulge for —10° < / < 10° and 
b = ±1° (GMV12) as weU as for b = -5°, where the 
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Fig. 1. — Mean radial velocity and velocity dispersion versus lon- 
gitude for the model's boxy bulge, as seen from the Sun for different 
latitudes |6| = 4°, 6°, 8°. Note how the dispersion decreases with 
latitude, and that in this projection the rotation is only approx- 
imately cylindrical; i.e., the slope of the velocity curve increases 
slightly towards the plane, as also seen in the Galactic bulge data. 

predictions of t h e mod el were subsequently confirmed by 
iGonzalez et alJ (|2012| ). Here we show that the model's 
boxy bulge at that time also shows approximate cylin- 
drical rotation, with the slope of the mean velocity curve 
along / increasing slightly towards lower latitudes; see 
Fig. [H This is sim ilar to the BRAVA kine matic data 
( Howard et alJ 120091 ) and the simulation by iShen et all 
(|2010( ) . The cylindrical rotation in boxy bulges is clearest 
in edge-on view, as shown in previous simulations (e .g., 
iCombes et "allll990t lAthanassoula fc MisiriotisI 12002'): it 
is weakened by the nearly end-on orientation of the Milky 
Way bar and the perspective effects. As in our previous 
work (MVGll, GMV12), we here take the solar galac- 
tocentric distance Dq = 8 kpc, scale the bar length to 
4.5 kpc, and use a bar orientation of a = 25° with respect 
to the Sun-Galactic center line. 

We study the formation of bulge metallicity gradients 
in the unstable disk - boxy bulge scenario with a simple 
illustrative model. We assume that the initial exponen- 
tial disk follows a specified radial metallicity gradient, 
which we imagine is set up during disk build-up prior 
to bar formation. We assign to each part icle a metal- 
licity depending on its initial radius (e.g., iFriedli et al.l 
[l99l . according to [M/H]{R) = [M/H]o + aR x i?/kpc. 
After some explorati on we choose \M/H]n = 0.6 and 
an = -0.4. Unlike iBekki fc Tsuiimotol (|2011| ). we do 
not link these metallicities to the present-day Galactic 
disk near the Sun, because the buckling instability in 
the Milky Way must have happened long ago when the 
outer disk may have been only incompletely assembled. 

3. RESULTS 

3.1. Vertical metallicity gradient in the Milky Way bulge 

The bar and buckling instabilities leading to the for- 
mation of the boxy bulge strongly change the particle or- 
bits. Thereby the initial disk metallicity distribution is 
mapped into the bulge accordingly. Fig.[2]shows longitu- 
dinal and vertical metallicity profiles for the bulge at Tf,. 
The vertical gradient for 3° < \b\ < 7° along the minor 
axis is ab = —0.064 dex/deg, compared to the observed 
slope in the Milky Way bulge aMW = —0.075 dex/deg 
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Fig. 2. — Longitudinal and vertical metallicity profiles of the 
model as seen from the Sun. Left: with I for \b\ = 2°, 4°, 8°; 
the gradient along I decreases towards larger distances from the 
Galactic plane. Right: with |6| for I = 0°,5°,10°; the vertical 
gradient decreases with distance from the minor axis. 

from Zoccali et al." (20081) and ab = —0.06 dex/deg 
from Gonzalez ct al. (2 01ll ). Near the Galactic plane, 
|6| < 3°, the vertical gradient becomes flatter and even 
positive. The positive gradient is due to the contami- 
nation from foreground/background disk particles. The 
flat part comes from the mixing during the buckling in- 
stability in the inner regions that evolve into the cen- 
tral near-spheroidal component (GMV12), and also from 
the greater fraction of low-latitude outer bar particles. 
The absence of a clear vertical gradient near the Galac- 
tic plane is c onsistent with recent results from lRich et al.l 
()2007l [20121 ). However, their measured mean metallici- 
ties for M-giants ([Fe/H]~ —0.2) are slight ly lower in 
Baade's window than those for K-giants fjZoccali et al.l 
l2008[) and for our simple model. 

Fig. [3] shows the resulting MDFs in several minor axis 
flelds in the boxy bulge at Tf, = 1.9 Gyr. Note the clear 
gradient of mean metallicity with latitude, and that sim- 
ilarly to the observed histograms, it arises mostly from 
the decreasing fraction of higher metallicity stars away 
from the Galactic plane. The maximum metallicity in all 
histograms is ~ 0.6, the central value in the initial disk, 
and the lowest value is ~ —1.2. 2 Gyr later, at T;, = 4 
Gyr, the model metallicity gradient is nearly identical, 
while the mean metallicities in the minor axis fields have 
slightly decreased, due to the capture of lower- metallicity 
disk stars by the slowly growing bar. Had we assumed a 
shallower radial gradient in the initial disk, also the final 
vertical gradient in the bulge would be smaller. 

3.2. Full metallicity map for the boxy bulge 

So far we have shown that a radial metallicity gra- 
dient in the unstable initial disk can survive through 
the bar and buckling instabilities, and generate, for suit- 
able parameters, a vertical metallicity gradient similar 
to that observed in the Milky Way bulge. For compari- 
son with upcoming survey results (e.g., VVV, ARGOS, 
Gaia-ESO), we now provide a full metallicity map, which 
gives a different view of the predicted metallicities for 
this scenario. Differently from Fig. [H we exclude fore- 
ground/background disk stars with a distance cut, using 
only particles with 4 < _D < 12 kpc from the position of 
the Sun. The assumed initial disk metallicity distribu- 
tion and simulated galaxy snapshot are as before. Fig. 2] 
shows the resulting average bulge metallicity map over 
the area on the sky extended by the Galactic bulge. It 
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Fig. 3. — Metallicity histograms (MDFs) for model particles in 4 fields along the minor axis of the boxy bulge {left panel s). The mean 
value for each MDF is given by the red arrow pointing up from the bottom. The corresponding mean value for the data oflZoccali et ahl 
((2008) is indicated by the blu e arrow from t he top. The histograms in the right panels are based on data from [Zoccali et an ~ if2008l ) for 
b = -4°, -6°, -12°, and from lJohnson et al.l ||2011,) for b = -8°. 



has several noteworthy properties: 

(i) The approximate outline of the boxy bulge can 
be seen together with low-metallicity indentations in 
the Galactic plane. The latter are due to fore- 
groud/background stars in the planar bar, which in these 
positions dominate the bulge stars. 

(ii) A metallicity gradient is present in all directions 
on the sky, both vertically and radially. This is expected 
from the binding energy argument of Sect. 13.31 How- 
ever, in the central few degrees the gradients become 
shallower. 

(iii) Iso-metallicity contours are more elongated verti- 
cally than horizontally, whereas the surface density con- 
tours are flattened to the plane (MVGll). 

(iv) The asymmetry between I > and Z < at inter- 
mediate latitudes is due to perspective effects and signi- 
fies the bar origin. In Fig. |4] high metallicities extend to 
larger \l\ on the I > side at |6| ~ 5°, but for |6| ~ 10° 
they extend to larger |Z| on the I < side. 

This map has some remarkable similarities with that 
being constructed fro m the near-infrared photometric 
VISTA VVV survey (jGonzalez et al.l I2013D . For now 
we can compare quantitatively with a few more mea- 
sured values in off-axis fields. Details will of course de- 
pend on the parametrisation of our simple model. In 
the longitudinal direction the gra dient along b = 2° is 
ai = -0.05 dex/deg (see Fig. ^. iMinniti et all ()1995f ) 
found mean metallicities [Fe/II]= —0.3 in two fields at 
{l,b) = (8°, -7°) and (10°, -7°) for bulge stars selected 
with [Fe/H]> —1. The model values in these fields are 
-0.30 and -0.28. 
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Fig. 4. — Metallicity map of the model bulge and bar in galactic 
coordinates. Foreground and background disk particles with dis- 
tances < 4 kpc and > 12 kpc from the solar position are excluded. 
The colour in each square corresponds to the average metallicity 
in a cone with radius 0.5° centered at this position. 



3.3. Origin of the vertical gradient 

Clearly, while the bar and buckling instabilities scram- 
ble the orbits of disk stars, they do not do so enough 
to completely erase the preexisting metallicity gradient. 
High-metallicity stars tightly bound to the Galactic cen- 
ter initially remain more tightly bound in the final bulge, 
and initially more loosely bound stars with lower metal- 
licities end up at larger final radii, on average. This can 



be quantified by considering the change in Jacobi energy 
Ej^^v^ + ^R,cf>,z)-^nlR^ (1) 

in the rotating frame of the boxy bulge and bar during 
the evolution. Here we use standard cylindrical coor- 
dinates, V is total velocity, $ is the gravitational po- 
tential, and ilp is the pattern speed at Tj, — 1.9 Gyr. 
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Fig. 5. — Change of Jacobi energy between different times during 
the evolution of the model, evaluated in a rotating frame for pattern 
speed f!p(Ti,) = 40 km/s/kpc. For each bin in initial Jacobi energy 
i?jO , the mean value of Jacobi energy at time t is plotted with error 
bars denoting the standard deviation of the distribution of Ej{t) 
in the bin. Black points: EjQ(t = 0) versus Ejq (initial scatter 
in each bin). Blue points: -Ejfinal(* = TJ)) versus Ejq (final boxy 
bulge compared to initial disk). Red open circles (slightly displaced 
horizontally for clarity): Ej\^^^(t = Tbar) versus Ejo (after full bar 
growth but prior to buckling instability, compared to initial disk). 



Fig. [S] shows that particles are scattered in Jacobi en- 
ergy by the two instabihties, but only over a fraction 
of the available total range in I.e., they retain 

some memory of their initial values. It has also been 
shown that most bar particles inside the vertical in- 
ner Lindblad resonance (VILR) mix during the buck- 
ling instability, but stay in the inner bar regions, while 
most bar particles around VI LR are scattered to orbits 
which can visit larger heights ([Pfenniger fc Friedlilll991l : 
iMartinez- Valpuesta, fc Shlosmanll2Q04D . 

From these arguments follows that the number of high- 
metallicity particles in the final boxy bulge will decrease 
with height above the plane where i5j becomes less neg- 
ative. In the histograms of Fig. [U therefore, the number 
of stars on the metal-rich part of the distribution de- 
creases with latitude. The maximum metallicity in all 
histograms, but that at 6 = —12°, is still given by the 
assumed maximum metallicity at the center of the initial 
disk, here -1-0.6, because a small fraction of the tightly 
bound particles is scattered up even to 6 = 8°. 

The low-metallicity tail, on the other hand, is more 
prominent at high latitudes, because of the larger frac- 
tion of particles coming from the outer parts of the bar. 
Note that within the assumed model, a lower limit to the 
metallicity in the bulge fields is expected, which is set by 
the outer boundary of the part of the disk which partic- 
ipates in the instability. From the histograms in Fig. [31 
this is ~ —1.2, which is due to particles in the initial disk 
at i? = 4.5kpc, just inside the corotation radius of the 
bar before buckling, 4.7 kpc. 



These properties of the model histograms are surpris- 
ingly similar to those of the Galactic bulge MDFs dis- 
cussed in the literature. If our scenario were the full 
explanation for the metallicity structure of the Galactic 
bulge, the observed range in the metallicity histograms 
could be used directly to estimate the central metallicity 
and the radial gradient in the inner Galactic disk prior 
to bar formation and buckling. The observational his- 
tograms, however, also contain metal-poor stellar halo 
and thick disk stars, and possibly metal-rich disk stars 
formed after the buckling, which must be taken into ac- 
count in this argument. Finally we note that a radial 
metallicity gradient in the precursor disk could have been 
accompanied by a stellar age gradient. In the context of 
the present model, depending on the time of buckling, 
this could lead to a measurable vertical age gradient in 
the boxy bulge which would provide an independent test 
of the model. 

3.4. Boxy bulges in other galaxies 

Vertical gradients have so far been published 
only for a small number o f gala xi es with boxy 
bulges. iFalcon-Barroso et al.l (|2004l ): iJablonka et al.l 
(20071) found a vertical metallicity gradie nt in NGC 7332 
comparable to that in the Milky Way. iWilliams et al.l 
(2011) determined vertical gradients in two other boxy 
bulges, NGC 1381 and NGC 3390. They integrated along 
slits parallel to the major axis to increase signal-to-noise, 
and obtained a single bulge measurement at several z. 
A similar averaging over I G [— 5°,5°] in our model re- 
duces the measured gradient from az — —0.46 dex/kpc 
to az = —0.33 dex/kpc, because radial and vertical 
gradients are comparable. Still, in NGC 1381 the ob- 
servations show a strong gradient of averaged — 
—0.27 dex/kpc, while only a weak gradient is seen in 
NGC 3390, az = -0.1 dex/kp c (z e [5,10] arcsec). 
iPerez fc Sanchez-Blazaue3 ()2011[ ) found negative metal- 
licity gradients along the major axis for the maj o rity o f 
bulges in barred galaxies, while IWilliams et al.l (|2012D 
found a range of major axis metallicity gradients in boxy 
bulges, from negative to positive, with shallower sam- 
ple average than for unbarred early- type galaxies. In 
summary, the limited data suggest a range of metallicity 
gradients in boxy bulges, with the metallicity gradient 
in the Milky Way on the steep side. In the framework 
discussed here, this suggests a range of initial disk metal- 
licity gradients before buckling in these galaxies, perhaps 
depending on the time of bulge formation. 

3.5. Comparison with previous simulations 

The results presented in the previous subsecti ons are 
consistent with early N-body simulation work by ' Friedlil 
11993). He found that pre-existing vertical abundance 
gradients in the disk were quickly flattened both in the 
bar and disk regions but not entirely erased. Their 
models also included a shallow, outward radial gradient, 
which became much flatter in the disk region due to the 
mixing induced by the b ar, but was ap proximately pre- 
served in the bar region. iFriedlil (|1998l) also stated that 
when no initial vertical gradient was present, a negative 
vertical gradient appeared at i? — 0. However, in his 
model the initial gradient was an ~ —0.1 dex/kpc, while 
the model analyzed here has initial an — —0.4 dex/kpc. 



Vertical metallicity gradients 
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for similar bar size. In our model the final gradient along 
the bar is au = —0.26 dex/kpc and the final vertical gra- 
dient is = -0.46 dex/k pc (for b G [-3°, -7°]). Con- 
sistent with iFriedlil (|1998[ ) we also find that at larger / 
the vertical gradient is shallower (Fig[5J right). 

Recent N-body nu merical simulations by 
IBekki fc Tsuiimotol (|2QTTI) for a pure exponential 
disk initial model (PDS in their nomenclature) showed 
only a very shallow vertical gradient in the final boxy 
bulge. The difference to our results can be traced to 
the different initial disk metallicity profile. While their 
inner disk parameters are similar to ours, due to linking 
the profile to the present metallicity near the Sun, their 
profile becomes quite shallow beyond ~ 0.25 times the 
eventual final bar length. Thus the stars in the upper 
bulge which come from the outer parts of the bar are 
quite metal-rich. The difference, and the reason for the 
steep vertical gradient in the bulge of our model is its 
steep initial radial gradient in the disk. 

Finally we note that if the secular evolution proceeds 
slowly, a radial gradient in the disk before the instability 
may be (partially) erased by migration processes. Thus 
a steep final bulge gradient is favoured by rapid secular 
evolution such as in the model investigated here. 

4. CONCLUSIONS 

We can summarize our conclusions as follows: 

(i) The vertical metallicity gradient observed in the 
Galactic bulge can be reproduced with a secularly 
evolved barred galaxy model in which a boxy bulge 
formed after a buckling instability. In itself, a vertical 
gradient is therefore not a strong argument for the exis- 
tence of a classical bulge in the Milky Way. 

(ii) Mixing during the bar and buckling instabilities is 
incomplete, and therefore radial metallicity gradients in 
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the initial disk can transform into gradients in the boxy 
bulge. 

(iii) In this framework, the range of bulge star metal- 
licities at various latitudes constrains the radial gradient 
in the precursor disk. 

(iv) The full bulge metallicity map shows an overall ra- 
dial gradient on the sky as well as longitudinal perspec- 
tive asymmetries. Iso-metallicity contours are elongated 
vertically. 

(v) Depending on when the bulge formed and on the 
properties of the precursor disk at that time, boxy bulges 
may or may not show metallicity gradients. 

The success of our simple model in explaining Milky 
Way observations suggests that its underlying idea has 
some merit and should be pursued further. Future work 
also needs to consider chemical evolution and possible 
vertical metallicity gradients in the initial disk, star for- 
mation induced by the bar, the possible late build-up of a 
metal rich inner disk, a nd the contribution of halo stars 
and other components (jNess et al.ll2012bl ) to the bulge 
MDFs. 

Another important next step for understanding the ori- 
gin of the Milky Way bulge is to investigate the correla- 
tions between metallicities and kinematics in this model 
and in gener alisations of it. Comparison with similar 
observations (jBabusiaux et al.l l2010t iNess et all l2012bl ) 
is likely to shed light on the possible multi-component 
nature of the Galactic bulge. The results will be impor- 
tant also for interpreting the data for bulges in external 
galaxies. 
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